AD-A202  747 


_ UNCLASSIFIED _ f_  .  i» 

StCumTv  Cl»S4'»iC*tio>4  O'  f  "kOt  f*yi»n  D««  Ent*i*d} 


^  rno* 


REPORT  documentation  PAGE 


HEPonT  NUMBCn 

AIM  1078 


READ  INSTRUCTIONS 

_ I _ BEFORE  COMPLETING  FORM 

I.  OOVT  accession  NO.I  1.  RECIPIENT’S  CATALOG  NUMBER 


Title  f«nd  Subiiii9i 

An  Optimal  Scale  for  Edge  detection 


S.  TYPE  OF  REPORT  •  PERIOD  COVERED 

memorandum 


S.  PERFORMING  ORC.  REPORT  NUMBER 


Au  THORrt; 

Davi  Geiger  and  Tomaso  Poggio 


t.  contract  or  grant  numberi*; 

N00014-85-K-0214 


PERFORMING  organization  NAME  AND  ADDRESS 

Artificial  Intelligence  Laboratory 
545  Technology  Square 

Cambridge,  MA  02139 _ _ _ 

controlling  office  name  ano  aooress 

Advanced  Research  Projects  Agency 
1400  Wilson  Blvd. 

_ Arlington,  VA  22209 _ 

IF.  MONITORING  AGENCY  NAME  S  AOORCSS<l<  dlllmtanl  Inm  CanmlUng  Olll€9) 

Office  of  Naval  Research 
Information  Systems 
Arlington,  VA  22217 

IS.  distribution  statement  r«i  iRi*  R*pwi) 

Distribution  is  unlimited 


17.  distribution  statement  (•<  (N*  atsiracl  mfimd  In  BlpcB  SB.  If  Blffaranl  i 


10.  program  element.  PROJECT,  TASK 

AREA  A  WORK  UNIT  NUMBERS 


IS.  REPORT  DATE 

September  1988 

IS.  NUMBER  OF  PAGES 

^_23 _ 

rnrTECURiTYnELASi~7priMj"r;wIriir 


UNCLASSIFIED 


Is*.  DECLASSIFICATION/ downgrading 
SCHEDULE 


Unlimited 


DTIC 

,..;lecte| 

DEC  0  6  19881 


IB.  SUPPLEMENTARY  NOTES 


IS.  KEY  WORDS  (Caflilnua  •«  r«*a,a*  aids  If  nacaaaarjr  and  Idanllfr  Bp  Black  naatBarf 

edge  detection 

cross-validation 

regularization 


I ZO.  abstract  fCanilmia  an  ravaraa  alda  If  nacaaaarp  and  Idwnllfr  Bf  Black  n«aiBor,| 


See  back  of  page 


DD  1473  EDITION  OF  I  NOV  SS  IS  OBSOLETE 

s/N  o:es-oi4-<«oi  i 


UNCLASSIFIED 

security  CLASSIFICATION  OF  THIS  PAGE  fPAan  Oaia  Bnlara, 


8  8  12  5  0  7  2 


1  '  ‘  ^ 


Block  20  (cont'd) 


Many  problems  in  early  vision  are  ill  posed^.  Edge  detection  is  a  typical 
example.  This  paper  applies  regularization  techniques  to  the  problem  of  edge 
detection.  We  derive  an  optimal  filter  for  edge  detection  with  a  size  controlled 
by  the  regularization  parameter  A  and  compare  it  to  the  Gaussian  filter.  A 
formula  relating  the  signal- to- noise  ratio  to  the  parameter  A  is  derived  from 
regularization  analysis  for  the  case  of  small  values  of  A.  We  also  discuss 
the  method  of  Generalized  Cross  Validation  for  obtaining  the  optimal  filter 
scale.  FinaJly,  we  use  our  framework  to  explain  two  perceptual  phenomena: 
coarsely  quantized  images  becoming  recognizable  by  either  blurring  or  adding 
noise. 
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Introduction 


Edge  detection  can  be  considered  as  a  problem  of  numerical  differentiation 
which  can  be  shown  to  be  an  ill-posed  problem  One  approach  to  solve  such 
an  ill-posed  problem  is  to  regularize  it.  Standard  regularization  techniques 
suggest  the  use  of  Gaussian-like  filters  before  differentiation^’'*.  In  this  paper, 
we  address  the  issue  of  how  to  estimate  the  optimal  scale  of  the  filter,  that  is, 
the  amount  of  smoothing  required  by  the  given  image  data.  More  precisely, 
we  will: 

1.  Obtain  a  regularization  filter  for  edge  detection  and  show  that  it  can 
be  seen  as  a  generalization  of  the  Poggio,  Voorhees  and  Yuille  filter 
(1985)  (PVY  filter)  ^ 

2.  Understand  the  role  of  the  parameter  A  (a  problem  mentioned  by 
T. Poggio  and  V.Torre  1984  ^)  and  find  out  how  to  calculate  it.  We 
will  show  that  the  optimal  value  of  A,  which  controls  the  size  of  the 
filter,  is  related  to  the  signal-to-noise  ratio. 

3.  Compare  the  optimal  filter  with  the  Gaussian  filter  and  PVY  filter.  We 
show  how  stable  <t  (the  size  of  the  gaussian  filter)  is  with  respect  to  A. 
This  also  corresponds  to  the  stability  of  the  width  of  the  optimal  filter 
with  respect  to  the  parameter  A. 

4.  Show  that  the  perceptual  phenomena  of  improved  recognizability  of 
coarsely  quantized  images  by  either  adding  noise  or  blurring  can  be 
understood  in  terms  of  scale  changes  in  the  edge  detector’s  filter. 

5.  Apply  this  filter  to  edge  detection.  Given  the  image  data,  we  use 
Canny’s^  edge  detector  but  using  the  scale  provided  by  the  regulariza¬ 
tion  analysis.  More  precisely,  the  noise  in  the  image  is  first  computed, 
then  the  signal-to-noise  ratio  is  obtained  and  finally  the  A  parameter 
is  used  to  set  the  scale  of  the  filter  used  in  computing  the  edges. 

6.  Compare  the  standard  regularization  techniques  with  the  Generalized 
Cross  Validation  method  for  finding  the  optimal  A. 

7.  Propose  biological  mechanisms  for  selecting  the  optimal  scale. 
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We  now  give  an  outline  of  the  paper.  Chapter  2  is  a  brief  summary  of 
results  relating  to  ill-posed  problems,  regularization  analysis  and  the  varia¬ 
tional  method.  In  Chapter  3  we  derive  the  optimal  filter.  We  show  that  this 
filter  gives  the  solution  of  the  ill-posed  problem  of  finding  the  exact  signal 
from  noisy  data.  In  Chapter  4  we  compare  the  optimal  filter,  Gaussian  filter 
and  PVY  filter.  In  Chapter  5  we  analyze  the  role  of  the  parameter  A.  In 
Chapter  6  we  present  an  implementation  of  the  edge  detection  results  from 
Chapter  4  for  computational  vision.  In  Chapter  7  we  give  a  possible  expla¬ 
nation  for  the  phenomena  of  improved  recognizability  of  coarsely  quantized 
image  when  either  noise  is  added  or  blurring.  In  Chapter  8  we  compare 
the  method  of  Generalized  Cross  Validation  (GCV)  for  finding  A  with  the 
standard  regularization  analysis. 


2  Framework  for  Edge  Detection 


Regularization  Techniques  for  Hi-posed  Problems 
X  problem 

Az  =  u  (2.1) 

for  which  the  class  5  of  solutions  z,  given  A  and  u,  is  not  compact  (changes 
on  the  right-hand  side  of  the  equation  can  take  u  outside  the  set  AS)  is 
called  ill-posed.  The  approach  suggested  by  Tikhonov  to  deal  with  ill-posed 
problems  is  to  construct  approximate  solutions  of  equation  (2.1)  that  are 
stable  under  small  changes  in  the  data  u. 

If  the  right-hand  side  of  equation  (2.1)  is  known  only  approximately,  we 
have  u{x,y)  =  urix^y)  +  v{x,y)  ,  where  UT{x,y)  is  the  true  solution  and 
v{x,y)  is  noise.  Then, 


where 


z(u;,i/) 


fc(u>,i/) 


zt{ui,  v) 


v{u,v) 

k{u},u) 


Az  =  [  f  K{x  -  i,y -T)z{i,T)didT  (2.2) 

J  —  OO  J  — oo 

and  is  the  Fourier  transform  of  K{x.,y).  It  would  seem  natural  to 

take  the  solution  of  equation  (2.1)  as  being 
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since  ux{u},iy)  =  k{uj,i>)zx{u},i').  However,  this  function  may  not  exist  since 
the  last  integral  may  diverge.  Furthermore,  even  if  this  ratio  does  have  ?n 
inverse  Fourier  transform,  the  deviation  from  zero  (in  the  C-  or  £2-nietric) 
can  be  arbitrarily  large,  and  thus,  we  cannot  think  of  the  exact  solution  of 
equation  (2.1)  as  an  approximate  solution  of  the  equation  with  approximate 
right-hand  side. 

Finding  edges  in  an  image  is  in  general  an  ill-posed  problem'*,  since  it 
involves  taking  an  appropriate  derivative  of  noisy  data  (notice  that  we  do 
not  specify  which  derivative  operator  should  be  used:  it  may  be  a  directional 
derivative  *  or  any  other  desirable  differential  operator).  The  differentiation 
of  the  function  u(i)  is  ill-posed,  since  it  can  be  seen  as  a  solution  of  equation 
2.1  for  the  operator  A  of  the  form 

/  h{x  -  i)z{S,)di  =  u{x) 

J -OO  J -oo 

where  h{x)is  the  step  function.  As  described  by  Rheinsch^  and  by  Poggio  and 
Torre  this  problem  can  be  regularized  by  smoothing  the  data  before  taking 
derivatives.  The  idea  is  to  consider  the  regularized  solution  z  of  equation 
(2.1),  with  A  being  the  imaging  operator,  such  that  z  is  sufficiently  well- 
behaved  for  numericed  differentiation  and  as  close  as  possible  to  the  true 
data. 

Tikhonov^  proves  that,  for  the  case  of  one  dimensional  image  data,  to 
approximate  a  solution  of  equation  (2.1)  one  takes  the  solution  of  a  different 
problem,  the  one  of  minimizing  the  functional  given  by  the  following  equa¬ 
tion,  that  is  close  to  the  original  problem  for  small  values  of  the  error  in  the 
data: 


M^[z,u]  =  f  f  [Az  —  dx  dy  \tl[z\  (2.3) 

J  —  OP  J  —  OO 

where  u{x,y)  is  the  image  data,  A  is  the  regularization  parameter,  A  is  in 
our  special  case  a  convolution  operator  and  fl[z]  is  the  stabilizing  operator. 
We  will  be  considering  the  special  case  where 
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f<x>  foo  /  \ 


dxdy. 


The  operator  fl[2|  forces  to  smooth  the  data  and  the  motivation  to  choose 
this  particular  stabilizing  operator  is  to  guarantee  a  desired  smoothness  in 
the  solution.  The  order  of  the  derivative  in  nfzj,  controlled  by  the  parameter 
p,  should  be  high  enough  to  ensure  the  appropriate  degree  of  differentiability 
in  z  required  by  the  derivative  operations  that  has  to  be  performed  next. 

It  should  be  pointed  out  that,  whereas  the  original  problem  (2.1)  does  not 
have  the  property  of  stability,  the  problem  of  minimizing  the  functional 
u]  is  stable  under  small  changes  in  the  right-hand  side  u.  This  stability 
has  been  attained  by  narrowing  the  class  of  possible  solutions  by  introducing 
the  stabilizing  functional  fifz]. 


3  The  Optimal  Filter 

The  Fourier  transform  of  (^  +  is  M(u;,  i/)2(u;,i/),  where  = 

(u;2  -Jr  u^Y'  Using  Parseval’s  theorem  we  can  rewrite  equation  (2.3)  as 


f  fOO  fOO 

-f- A  /  f  (w* -f  r/)2(— u»,  — 

J  —  OO  J -oo  J 

The  associated  Euler- Lagrange  equation  is 


dM^ 

dz{—uj  —  J/) 


[fc(u;,  i^)z(ui,  u)  —  u{u,  t/)]k{—u;,  —i/)  -f  A(u;*  +  i/*)'’2(u;,  i/)  =  0. 


For  the  special  case  where  p  =  2  (stabilize  with  Laplacian)  and  the  operator 
A  is  given  by  (2.2)  with  k{u>,u)  =  e“ we  obtain  the  regularized 
solution 


1  fOO  fOO 


oo  1  -b  A(u;*  -b 


2(u;,  i/)e-‘<“*+‘'>'>du;di/ 


,  — i/)]dci; 
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Optimal  Filter 


Figure  1:  The  optimal  filter  is  plotted  in  the  space  domain  for  different  values 
of  the  regularization  and  scale  parameter  \ 

that  corresponds,  in  the  Fourier  domain,  to  filtering  the  data  with 

/(u;,i/,A)  =  ^  (3.1) 

The  term  k{u;,i/)  approximates  the  point  spread  function  of  the  imaging 
device.  For  the  human  eye  the  fc(u>,  u)  is  described  by  a  Bessel  function  and 
a  Gaussian  can  be  a  good  approximation.  In  this  case  the  bandwidth  (<t) 
of  the  point  spread  function  of  the  pupil  is  of  the  order  of  20  seconds  of  arc 
which  implies  b  =  0.5<t^  =  1.5  10“®  per  degree  square.  For  CCD  cameras 
the  bandwidth  is  less  than  a  pixel,  restricting  the  effect  of  this  term.  In  our 
experiments  we  used  b  =  0.5  per  pixel  square. 

The  case  b=0.0  gives  the  PVY  filter.  In  this  case  the  filter  is  not  smooth 
enough  to  ensure  differentiability  with  a  second  order  differential  operator 
such  as  the  Laplacian  and  a  higher  order  stabilizer  is  needed  For  6  >  0 
this  problem,  however,  disappears.  The  filter  is  plotted  (6  =  0.5)  in  the  space 
domain  for  different  values  of  the  parameter  A  (see  figure  1),  which  controls 
the  size  of  the  filter. 

Appendix  2  shows  how  to  obtain  the  Wiener  filter  using  the  variational 
method.  We  can  see  that  the  Wiener  filter  is  slightly  different  from  the 
optimal  filter,  because  of  the  stabilizer  term. 
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4  Filters:  Optimal,  Gaussian,  PVY  and  Wiener 

4.1  The  relation  between  A  ,  a  ,  a  and  7 

In  order  to  compare  the  filters,  we  have  to  assume  an  equivalent  criterion  for 
the  Optimal  filter,  the  Gaussian  filter,  PVY  filter  and  Wiener  filter.  A  pos¬ 
sible  criterion  is  to  choose  the  same  half  amplitude  in  the  frequency  domain. 

1)  Gaussian-filter  -  1  then  u;-  = 

2)  Optimal  filter  -  j  ,  then  Aw*  =  e"*""* 

3)  PVY-filter  -  \  then  aw*  =  1 

4)  Wiener- filter  -  A  =  then  7  = 

Then  we  plot  a  graph  <t  x  A  for  b=0.5  per  pixel  square  (see  figure  2) 

We  notice  here  the  stability  of  the  parameter  a  (the  bandwidth, in  pixel 
unity, of  a  filter)  with  respect  to  the  parameter  A. 

Another  criterion  for  equivalence  is  to  chose  filters  that  have  the  same 
bandwidth  in  their  second-derivative  (see  figure  4).  This  is  important  when 
computing  zero-crossings,  since  in  this  case  the  desired  operator  of  convolu- 
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L  ^ 


tion  is  the  second  derivative  of  the  filter.  The  value  of  the  parameter  A  may 
vary  with  these  different  definitions,  but  the  actual  bandwidth  of  the  filter 
will  not  vary  significantly. 

5  The  Optimal  Scale 

In  order  to  find  the  optimal  scale  we  vary  the  parameter  A  to  obtain  the 
closest  solution  to  the  true  solution.  More  precisely  (see  appendix  3),  we 
minimize 


nri  /  X  /  M21  1  /•“  A^(u;*  + +  e  j  , 


where  £'{.}  is  the  expectation  value  operator,  S{u,i/)  is  the  spectral  density 
(power  spectrum)  of  zrix^y)  and  is  the  spectral  density  of  v{x,y), 

assuming  that  v(i,  y)  is  stationary. 

Differentiating  equation  5.1  with  respect  to  A  we  obtain 


Jo  Jo  +  + 


=rf 


*  N{u> 


Jo  Jo  (e-‘(“'+‘'*)  +  A(a;2  +  t/2)2)3  Jo  Jo  (e-‘-(“^+‘'*)  +  A(u;2  +  r/2)2)3 

(5.2) 

For  the  discrete  case  the  space  is  defined  in  a  grid  with  N  x  M  pix¬ 
els  and  the  spatial  frequency  given  by  Wi  =  and  t/j  =  where  i  = 

—  y,  ...,0, ...,  y  and  j  =  — y,...,0,...,  y.  We  follow  the  same  calculations 
as  in  the  continuous  case  and  keeping  the  order  of  the  stabilizer  a  general 
parameter  p  (before  we  set  p  =  2).  Using  the  symmetry  with  respect  to  zero 
of  5(u;,i/)  andN{ij}^u)  one  obtain  similarly  to  equation  5.2  the  equation 


du!  du 


i=oj=o  (e-“'((^)’+(A)’)  +  A(2jr)2i>((i)2  + 

N(i,j)e-‘'^-»*)’-^(AI’>(2ir)»((l)»  +  (ijYf 
i=o  j=o  ^e-bj,((A)>+(i)>)  +  A(2>r)2P((i)2  + 
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Within  the  range  of  the  parameter  A  there  are  two  cases  two  consider. 
One  for  very  small  values,  i.e.  when  A  — >  0,  and  the  other  case  for  A  not  small. 
The  second  case  is  applicable  whenever  the  signal-to-noise  ratio  is  small  and 
to  find  A  one  has  to  solve  equation  5.2b.  The  first  case  is  applicable  whenever 
the  signal-to-noise  ratio  is  very  small  and  the  following  short  cut  to  find  A 
can  be  used.  We  assume  the  asymptotic  values 


lim  S{jj,  u)  = 


>  0 


lim  ^V(u;,  u)  = 


No 


(u;^  -f-  u'^Y 


a  >  0 


and  considering  the  case  with  white  noise  (a  =  0)  we  obtain  (see  appendix 
4)  A  as  the  solution  of 


j\clp(^)7^J^26(^A)-Pr(4p-l))  /  (46(^A)c-p+2p+l)  . 

(5.3) 

For  the  case  b  =  0.0  (the  PVY  filter)  we  get  A  =  (2p  -f  l/4p  -  l)*^(^)c . 

A  graph  of  signal-to-noise  ratio  versus  A  for  6  =  0.5  per  pixel  square  and 
is  plotted  in  figure  3.  Thus  the  optimal  size  of  the  filter  can  be  obtained 
directly  from  an  estimate  of  the  signal-to-noise  ratio  for  any  given  image. 
The  equation  5.3,  for  finding  the  parameter  A,  turns  out  to  be  sensitive  when 
we  fix  the  value  of  the  signal-to-noise  ratio  and  vary  the  values  of  6  or  c. 
Hov/ever  the  relevant  parameter,  namely,  the  width  of  the  filter  is  not  very 
sensitive  to  it.  In  our  experiments  the  width  did  not  vary  much  from  0.8 
pixels  to  1.1  pixels.  The  parameter  c  was  varied  from  2  to  15,  the  parameter 
p  (degree  of  smoothness)  was  varied  from  2  to  6  and  the  parameter  b  (optical 
blurring)  was  varied  from  0  to  2.  We  have  to  remember  that  here  we  are  just 
considering  large  values  of  signal-to-noise  ratio,  typically  >  100. 


6  Noise  Estimation 

Here  we  assume  the  noise  to  be  caused  by  the  deficiency  in  the  optical  sys¬ 
tem  or  by  artificially  adding  a  random  variable  to  the  image.  So  in  this  way 
noise  is  independent  of  the  recognition  task  and  independent  of  the  scene  or 
a  particular  scale.  To  estimate  the  noise  we  use  a  technique  developed  by 
H.Voorhees  First,  the  gradient  of  the  image  is  computed.  A  Rayleigh 
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c 


=  10 


p  =  2 


c  =  3 


qna  1  -r»o  ’ 


P  =  3  c  -  10 

1  3  1  -'--j  1  5^ 

r* . .  ^ 


•  fi4  n  iv,  II  Oft  0  ft? 


Figure  3:  Signal-to-noise  ratio 
x-axis  is  given  by  the  width  of  ti 


p  =  2  c  =  JO 

5  t -jna  >  -  fio  I 
40ftft^~ 


JfiOft 


I)  -ift  ft  91  n  92  n  9  ? 

versus  for  b  =  0.5.  The  last  curve  the 
'.e  filter  cr 


distribution  is  then  fit  to  the  histogram  of  the  norm  of  the  gradient  and  the 
noise  parameters  are  estimated.  The  signal  power  is  obtained  from  the  vari¬ 
ance  of  the  intensity  of  the  image.  With  this  method  the  program  estimates 
from  the  image  data  the  corresponding  signal-to-noise  ratio.  This  ratio  gives 
A  from  relations  such  as  equation  5.3  or  5.2b.  The  results  are  shown  in 
figures:  4,  5  and  6.  One  alternative  for  noise  estimation  could  be  done  by 
taking  a  sequence  of  images  and  from  the  temporal  correlation  obtain  the 
noise  parameters.  Another  one  is,  assuming  the  noise  to  be  white,  to  find 
the  average  value  for  large  frequencies  of  the  image-spectrum,  since  typically 
for  large  frequencies  the  noise  is  constant  and  the  signal  values  close  to  zero. 


7  Two  perceptu€d  phenomena:  explanations 
and  biological  implications 

7.1  Coarse  quantized  in\ages  can  be  better  recog¬ 
nized  when  noise  is  added 

We  first  discuss  the  perceptu2d  phenomena  of  improved  recognizability  of 
coarse  quantized  images  when  noise  is  added^.  Consider  an  image  with  320 
by  384  pixels  and  8  bits.  A  coarsely  quantized  version  of  it  is  shown  in 
figure  5a.  The  optimal  filter  for  figure  5a,  estimated  as  explained  above. 
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Pom*®  gauMMn  2)optlmal-fitter 


Figure  4;  Zero- crossings  ,  the  Laplacian  applied  to  an  image  that  was 
smoothed  using  a  Gaussian  filter  and  the  optimal  filter  with  the  same  width 
for  the  second  derivative 

turns  out  to  have  a  small  scale  of  ^  1.0  pixel  corresponding  to  very  large 
signal-to-noise  ratio  of  So/ =  1000.0.  We  then  apply  the  Canny’s  edge 
detector  using  a  gaussian  filter  with  the  same  bandwidth  predicted  by  the 
regularization.  The  edges  do  not  easily  reveal  a  face  as  we  see  in  figure  5c. 
Gaussian  white  noise  with  standard  deviation  70  is  added  to  figure  5a  (see 
figure  5b),  making  recognition  easier.  Estimation  of  the  optimal  scale,  using 
equation  5.2b  since  the  signal-to-noise  ratio  is  small,  gives  now  a  width  of 
6.0  pixels.  The  corresponding  contours  (output  from  Canny’s  detector) 
reveal  the  face  in  a  much  better  way. 

These  results  may  shed  some  light  on  what  the  visual  system  may  be  doing. 
Harmon  and  Julesz  ®  claim  that  for  the  quantized  image  “high  frequencies 
introduced  by  quantized  blocking  mask  the  lower  spatial  frequencies  which 
convey  information  about  the  face,  preventing  recognition”.  In  our  frame¬ 
work  two  processes  determine  recognizability  of  the  face.  The  first  process 
consists  of  the  estimation  of  the  signal-to-noise  ratio  {So/ No).  The  second 
step  is  to  use  So/Ng  to  set  the  optimal  A  for  computing  the  corresponding 
“edges”.  In  the  case  of  the  quantized  image  the  ratio  So/ No  is  large.  A  is 
then  small,  which  implies  that  a  large  bandwidth  channel  (in  the  spatial- 
frequency  domain)  is  selected.  The  zero  crossings  for  this  channels  do  not 
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Figure  5:  a,c)  Quantized  image  and  edges  with  a  =  1.0  pixel.  b,d)  White 
noise,  standard  deviation  of  70  units,  was  added  to  the  quantized  image. 
Now  the  edges  are  computed  using  (T  =  6.0  pixels. 
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easily  allow  face  recognition  because  they  mostly  capture  the  box  outlines. 
For  noisy  quantized  images  the  ratio  Sg/ No  is  small  and  correspondingly  A  is 
large.  This  implies  a  filter  with  small  bandwidth.  In  this  case  the  small  band¬ 
width  filter  suppresses  the  noise  and,  as  a  side  effect,  also  the  high  frequency 
outlines  of  the  boxes. 

This  explanation  is  not  in  contrast  with  the  one  given  by  Canny  ^  or  by 
the  one  given  by  Morrone,  Burr  and  Ross  (1983)  when  they  claim  “that  added 
noise  (more  high  frequencies)  destroys  the  propensity  to  organize  the  image 
according  to  its  spurious  high-frequency  structure,...”,  but  is  more  precise. 


7.2  Improved  recognizability  of  coarse  quantized  im¬ 
ages  by  blurring 

Blurring  coarsely  quantized  images  also  improves  recognition*.  The  explana¬ 
tion  for  this  second  perceptual  phenomena  is  natural  in  our  model.  Consider 
the  zero  crossings  of  the  blurred  image.  For  the  Gaussian  filter  they  are 
represented  by  the  solution  of 

V^(G(<T2;i,t/)*(G(<Ti;i,y)*J(i,y)))  =  0 

Where  I{x,y)  is  the  quantized  image,  *  stands  for  convolution,  G{<T;x,y) 
is  the  Gaussian  and  the  blurring  process  is  {G{a’i;x,y)  *  I{x,y)).  Because 
all  of  those  operations  are  linear  we  can  first  convolve  both  Gaussians.  In 
the  Fourier  domain  this  is  equivalent  to  multiplying  the  Gaussians,  and  the 
multiplication  of  two  Gaussians  is  a  Gaussian  with  <r  =  [(tti  +  So 

computing  zero  crossing  of  a  blurred  image  is  the  same  as  computing  zero 
crossings  of  the  same  image  but  with  larger  (t  (<t  =  [((Ti  )’  +  (»,)>)■« ). 
Therefore  blurring  is  equivalent  to  using  an  effectively  larger  filter  for  edge 
detection.  This  has  the  effect  of  suppressing  the  spurious  high  frequency 
edges  introduced  by  coarse  quantization.  * 

^Notice  that  blurring  the  quantised  nouy  image  has  the  effect  of  increasing  the  esti¬ 
mated  signal-to-noise  ratio,  thereby  reducing  o-j  to  a  value  close  to  the  one  obtained  for 
the  quantized  image. 
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8  The  Method  of  Generalized  Cross  Valida 
tion  and  Regularization 


When  So/ No  cannot  be  directly  estimated,  it  is  natural  to  consider  the 
method  of  Generalized  Cross  Validation  (GCV)®.  The  GCV  method  states 
that  the  optimal  value  of  A  can  be  obtained  by  minimizing  the  functional 
(here  in  one  dimension) 


V(A)  = 


[AZn,x{ti)  -  U,]^ 

(l-afcfc(A))2 


(8.3) 


where  a;fc(A)  =  (1  -  afcfc(A))/(l  -  ^  Ej=o and  afcfe(A)  =  ^(Az„,A)(tfe). 
From  here  on  we  use  the  Gaussian  filter  ,  since  would  be  computationally 
more  expensive  to  use  the  optimal  filter.  Equation  8.3  reduces  to 


V(<t)  = 


n  (1 


2lt(T  ■' 


The  method  is  computationally  expensive  but  intrinsically  parallel.  We  have 
implemented  it  on  the  Connection  Machine  We  tested  this  method  on  dif¬ 
ferent  images  including  the  ones  in  figure  5  with  various  amounts  of  noise.  We 
notice  a  consistency  of  the  GCV  with  the  results  obtained  with  the  standard 
regularization.  The  GCV  gave  allways  (see  figure  6)  a  slightly  mailer  value  for 
the  parameter  A  (and  <r,  the  width  of  the  filter).  We  point  out  that  equation 
5.3  is  just  applicable  when  the  signal-to-noise  ratio  is  high  (higher  than  200) 
and  the  other  cases  we  used  equation  5.2b.  The  results  of  applying  Canny’s 
edges  detector  with  both  estimates  of  the  parameter  suggest  a  slightly  better 
performance  for  the  standard  regularization  method.  For  large  amount  of 
noise  as  in  figure  5b  we  obtained  significantly  different  results.  For  figure 
5c  for  example  the  GCV  method  gave  <r  =  3.5  pixels  as  oppose  to  6.0  given 
by  the  standard  regrilarization  method.  This  suggests  that  since  no  noise 
estimation  is  involved  in  the  GCV  method  whenever  the  signal-to-noise  ratio 
is  small  the  standard  regularization  is  going  to  give  better  results.  Typically 
slices  of  the  image  of  size  80  pixels  require  20  milliseconds  for  computing 
V(<r).  Using  Newton’s  method  to  find  the  minimum,  the  algorithm  con¬ 
verges  after  10  iterations  of  V{<r).  Therefore  the  GCV  method  takes  in  this 


^The  Connection  Machine  is  a  massivelly  parallel  computer  with  65,536  processors 
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case  0.2  seconds  to  find  the  optimal  a.  In  figure  6  we  show  several  examples 
where  the  optimal  scale  was  estimated  with  the  standard  regularization  and 

GCV'  methods. 

It  is  interesting  to  notice  that  for  the  figure  6. a  the  optimal  scale  was  for 
A  =  0.1  (<7  =  0.91  pixels),  however  there  are  two  distinct  scales  for  this  image 
that  could  be  considered  optimal.  The  coarser  scale  was  not  captured  with 
any  regularization  method.  A  possible  way  of  obtaining  both  scales  would 
be  by  reconsidering  the  definition  of  noise  or  signal.  For  instance  if  the  noise 
was  defined  to  be  trie  medium  and  high  frequency  of  the  spectrum  then  the 
noise  estimation  would  be  higher  than  the  one  we  used  and  the  coarser  scale 
would  be  optimal. 


Figure  6.  Sere  ml  eramples  of  applying  the  stan-  rrd  regularization  method 
o.nd  the  G'd  metixod  to  obtain  the  optimal  scale 


Westminst  and  Canny’s  edges  with  a  -  0.86  and  cr  ~  0.66.  Cross- 
validation:  (T  ^  O.titi  pizels,  standard-regularization:  <t  =  0.86  pixels 


C’loth  and  Canny’s  edges  with  a  =  0.91  and  a  —  0.54.  Cross- calidat ion 
<T  -  0.54  pixels,  standard-regularization:  a  =  0.91  pixels. 


Walter  and  Canny’s  edges  with  a  =  1.63  and  c  =  1.57.  Cross-validation: 
cr  =  1.57  pixels,  standard-regularization:  (T  =  1.62  pixels.  We  solved  equation 
5.b  since  So/ No  ratio  smaller  than  200. 


.Apple  and  Canny’s  edges  with  <t  =  1.2  and  <t  ~  0.66.  Cross-validation: 
(T  —  0.66  pixels,  standard-regularization:  cr  =  1.2  pixels.  We  solved  equation 
5.h  since  So/ No  ratio  smaller  than  200. 
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Ericqn  and  Canny’s  edges  with  a  =  6.00  and  <t  =  3.5.  Cross-validation: 

<7  =  3.48  pixels,  standard-regularization:  o  =  6.08  pixels.  We  solved  equation 

5.b  since  So/ No  ratio  smaller  than  200. 

9  Conclusion 

We  have  derived  rigorously  the  optimal  way  of  filtering  images  prior  to  nu¬ 
merical  differentiation.  We  also  obtained  the  precise  relation  between  the 
scale  of  the  filter  A  and  the  signal- to- noise  ratio  of  the  image.  Some  bio¬ 
logical  implications  were  also  considered.  In  particular  we  suggested  that 
humans  can  estimate  the  signal- to- noise  ratio  in  the  image  &om  which  the 
scale  A  is  computed.  Only  channels  with  the  appropriate  spatial-frequency 
band  are  then  used,  the  others  being  inhibited.  In  this  framework  it  is  possi¬ 
ble  to  understand  th'  perceptual  phenomena  of  improved  recognizability  of 
coarsely  quantized  images  when  noise  is  added.  When  the  signal-to-noise  ra¬ 
tio  is  large,  the  estimated  A  is  small  and  the  associated  edges  do  not  provide 
good  information  for  recognition.  When  So/ No  is  smaller,  the  estimated  A  is 
larger:  the  edges  then  provide  better  information  for  recognizing  the  face  in 
the  image.  When  the  signal-to-noise  ratio  cannot  be  estimated,  it  is  possible 
to  use  the  method  of  Cross  Validation  for  estimating  the  optimal  A. 
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D.G.).  Thanks  to  E.  Hildreth  for  the  suggestions,  to  E.  Gamble,  M.  Gen- 
nert,  and  H.  Voorhees  who  were  always  helpful  and  to  J.  Little  for  the  last 
comments. 


Appendix  1:  Deviation  of  the  regularized  solution  from  the  ex¬ 
act  one  and  the  stability  of  the  solution 

The  regularized  solution  can  be  written  in  the  form  (see  Chapter  3) 


/°°  v(u},  du 

(27r)2  y_oo  J-oo  K{lo)K{u)  ^  ’  ’ 


then 


z,(x,y)  -  zr(x,y)  = 


K{uj)Kiv) 


i—  f°°  ^  v(ij  du  (l.al) 

{2iry  J-oo  K{u)K{v)  '^  '  ^  ^ 


Now  one  defines 


As(<,A) 

A) 


1  /“  {/(g;,  t/,A)-l} 

(2ir)2  7-00  K{uj)I<{u) 

1  p  p  /(u;,t/,A) 

(25r)2  7-00  7-00  /C(u;)/<'(I/) 


(l.a2) 

(l.a3) 


so  that  z\{x,  y)  —  zt{x,  y)  =  ^s{^i  Vi  A)  +  A/^r(z,  y,  A)  One  notices  that  equa¬ 
tion  (l.a2)  characterizes  the  influence  of  the  regularization  and  equation 
(l.^l3)  characterizes  the  influence  of  the  noise. 

Under  these  assumptions,  the  variance  of  the  random  function  Ajv  is 
A)}  =  ^’(A)  =  jL/” 
and  A^(A)  =  suptAs{t,  A) 

How  stable  is  the  regularized  solution;  more  precisely,  how  much  does  the 
regularized  solution  change  with  a  change  in  the  order  of  the  stabilizer  p  ? 

It  is  possible  to  show  that: 

1)  The  asymptotic  (  u>,  i/  — »  oo  )  behavior  of  As  and  An  are  the  same  for 
stabilizers  of  the  form 


M{uj,i/)  =  {uj^+v"^Y+ap-i{u^+t/^Y~^+...+ao  or 


M(u;,v)  = 
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2)  From  formula  5.1 


T(A,p) 


2ir^\Jo  [L{ij},u)  +  u'^YY  ^•'0  [L{(jj,u)  +  u'^YY 


obviously,  A  is  a  function  of  p.  Consequently,  r(A,p)  =  V’(p)  and  the  function 
t/»(p)  has  a  unique  minimum  at  p  =  Po  =  c  —  a,  and  increases  monotonically 
in  the  interval  p  >  Po  with  a  finite  limit  for  p  — »  oo.  Tikhonov  ^  has  shown 
these  results  for  the  one  dimensional  case. 

Defining  the  measure  of  stability  for  arbitrary  p  >  0  as  being 


one  has  from  Tikhonov  ‘ 


Hv)  -  i’jPo) 
Hpo) 


0  <  ^(p)  ~  <  7o 

V'(Po)  1  -  To 

where  7<,  =  2n+2e-2a  ^  order  of  smoothness  of  the  operator  iir(u;) 

,  i.e.,  how  fast  k(<*>)  -*  0  when  (a  -*  oo  . 

For  the  case  where  k{tjj)  —  then  n  -+  oo.  Therefore  -*  0, 

which  is  to  say  that  the  function  ^(p)  is  weakly  dependent  on  the  order 
p  of  the  regularization.  Consequently,  one  can  quite  justifiably  replace  the 
optimal  order  of  regularization  po  with  an  arbitrary  order  p  >  Po.  For  white 
noise  Po  =  c  and  for  p  =  2  one  will  be  safe  whenever  working  with  images 
that  at  high  frequencies  do  not  fall  to  zero  faster  than 
It  is  clear  is  that  if  dealing  with  images  with  different  asymptotic  behavior, 
one  can  reset  p  to  give  a  good  solution,  and  the  stability  for  higher  p  will  be 
guaranteed. 

In  practice  we  are  interested  in  the  stability  of  the  scale  of  the  filter  and 
we  notice  that  the  actual  width  of  the  filter  did  not  vary  more  than  0.4  pixels 
when  the  parameters  c,  p  and  b  (the  point  spread  function  parameter)  were 
varied. 
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Appendix  2:  The  connection  between  the  regularization  method 
and  optimal  Wiener  filtering. 

The  regularized  solutions  of  equation  2.1  is 

(27r)^  J-oo  L{u),i/)  +  \m{u),v) 
and  can  be  written  in  the  form  of  a  convolution 

/OO  too  fOO 

/  /  Lx{t  -  T)u{T)dT 

■OO  j  — OO  */— OO 

where 


Lx{t  -  t) 


L{uf,i/)  +  AM(u>,  t/) 


Now  taking  T(XM)  (formula  5.1),  from  the  condition  that  this  functional  be 
minimized  on  the  set  of  functions  M(u?,  i/),  one  finds  by  elementary  calcula¬ 
tions  that  the  minimum  is  attained  with  the  function 


M  (w,  i/)  =  Mo(a>,  I/)  = 


X  u)  ‘ 


Therefore  the  approximate  (regularized)  solution  Zop(x,i/)  of  equation  2.1 
obtained  is 


^op' 


which  is  independent  of  the  parameter  A.  It  coincides  with  the  result  of 
applying  the  optimal  Wiener  filtering  to  find  z{x,  y)  from  «(x,  y)  =  ut{x,  y)  + 
v{x,y)  . 
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Appendix  3:  The  expression  for  E{\zx{t)  — 

We  show  the  result  for  a  more  general  class  of  solutions  where  the  smoothing 
term  is  given  by  A/ Notice  that  i/)  =  k{-ui)k[-v)k{ui)k{v). 

We  have 

zx{x,y)  -  ZT{x,y)  = 

(27rj  J—ooJ—oo  Z((w,  t'j  “h  AA/(w,  t') 

=ihr  r  ( = 

[ZTT J ~<x>  J —oo  L\Jjt}yU)  XM{u}^V) 

(2x)2  7-00  7-00 ■^I(u;,i/)  +  AM(a;,t/)  i(u;,  i/)  +  AM(u;,  i/)  i(u;®  + 

since  i/)  =  k{u)k{v)zT{u},u).  Therefore 
£^{\2\{x,y)  -  ZT{x,y)\^]  = 

I  (27r)*  7-00  v-oo  Lfw,*/)  +  AM(u;,  I/)  i 

■(2t)^J-^LJ  Kw-Xj  +  AMK-') - ^ 

“  (4x2)2  7-00  7-00  7-00  7-00 ^[L(u;,  I/)  +  AM(w,  i/)][L(u;',  I/O  +  AM(u/,  i/O] 

[L{u,  I/)  +  u')  +  i/')] 

since  E{v{u>,v)}  =  E{v{u)',i/')}  =  0.  For  stationary  random  processes, 

E{zt{u},  j/).zr(u;',  u')}  =  5(u;,  i/)6{u  +  u}')6{t/  +  i/') 

E{v(u,  i/)v(u;,  I/')}  =  N(u;,  i/)S(u;  +  u;')S(i/  +  i/'), 

where  S(u;  +  (*;')  is  Dirac’s  delta  function.  Performing  the  integration  in 
the  last  expression  with  respect  to  w'  and  using  the  properties  of  the  delta 
function  and  the  fact  that  M{w)  and  //(u;)  are  even  functions,  one  obtains 
the  value  of  the  deviation 

=  t(xm)  = 

4x2  [L{u),  I/)  +  XM{u;,  t/)]2 
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Appendix  4:  Solving  equation  5.2,  based  on  Tikhonov  ^ 

The  argument  for  using  the  asymptotic  expressions  (p  is  the  order  of  differ¬ 
entiation  in  the  smoothing  term)  follows  immediately  from 


no 

J 


A 


+  A(u;2 -f  ^ 


since  <p^(A)  — >  oo  as  A  — ►  0.  On  the  other  hand,  for  any  fixed  A  >  0, 
the  integral  <p^(A)  converges.  Therefore,  for  sufficiently  small  values  of  the 
parameter  A  ,  the  basic  contribution  to  tp  is  provided  by  large  frequencies  u. 
Then  one  can  replace  equation  5.2  by 


/■oo  yoo  ^  j  j  foo  roo  +  1/*)'’"“  j  j 

Jo  Jo  A(u;2  ~  Jo  Jo  +  + 

(5.26) 

with  the  change  of  variables  u)  =  ^  cos6  and  v  =  ^  sind  which  imply  = 

-I- and  duidu  =  ^  d^  dO  .  We  can  rewrite  equation  (5.2b)  as 


Jo  {Luo  +  H^r  ^~Jo  {LaM)  +  H^r  ^ 

where  Taj(0  =  The  integrals 


(5.36) 


^  roc  ^ 

"  Jo  {LUO  +  Jo  {LUO  +  HO' 

are  evaluated  by  the  method  of  steepest  descent  and  are  equal  to 


h 


.  [~jZZ 


+  0(At)) 


/,= 


+  Aff}’ 


2c+2a 


+  0(Ai)), 


where 


/i(^,A)  =  A/n 


{LUO  +  HO^ 
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and  the  double  prime  denotes  the  second  derivative  with  respect  to  Here, 
is  a  loot  of  the  equation 


{-2La,  +  -  (4p  +  2a)\r^L,.  +  =  0  . 

Substituting  the  values  found  for  Ii  and  I2  into  equation  (5.3b)  and  keeping 
only  the  principal  terms,  one  obtains  the  following  equation  for  determining 
the  optimal  A 


A  = 


-2p+2c-2a 
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